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SHOCK  WAVE  MODELLING  AND  THE  HULL  COMPUTER  CODE 


FORMATION  OF  MACH  STEMS 


1 .  INTRODUCTION 


Shock  waves  from  detonating  explosive  charges  are  capable  of 
inflicting  damage  over  wide  distances.  Their  formation,  their  propagation  and 
their  properties  at  any  instant  are  of  fundamental  importance  in  assessing 
their  effectiveness  against  various  targets.  Experimental  measurements  of 
shock  parameters  can  be  a  difficult,  costly  and  a  time  consuming  process.  Use 
of  appropriate  numerical  tools  is  probably  justified  when  estimates  are 
required;  confirmatory  experiments  can  then  be  reduced  to  an  absolute 
minimum. 


This  Report  describes  the  implementation  of  the  HULL  computer  code 
in  Australia,  and  its  initial  application  to  the  propagation  of  shock  waves  up 
an  inclined  ramp  in  a  Shock  Ttabe  and  the  study  of  any  associated  Mach  Stem 
formation  and  Triple  Point  trajectory. 


2.  HULL  CODE 


HULL  [1]  employs  an  Eulerian  finite-difference  scheme  which  can  be 
applied  to  either  two-or  three-dimensional  problems.  It  is  a  completely 
independent  package  made  up  of  three  main  programs  each  capable  of  running 
separately  plus  two  supplementary  programs.  These  programs  are:- 

KEEL  The  problem  is  defined  in  KEEL.  It  sets  up  the  calculation  grid  and 
stores  the  initial  values  for  each  cell  including  any  special 
boundary  conditions. 

HULL  The  number-crunching  is  carried  out  in  HULL.  Not  surprisingly,  this 
program  provides  the  name  of  the  whole  system.  Here,  the  files 
created  by  KEEL  are  used  in  the  finite  difference  scheme  to  advance 
the  problem  in  time  increments. 

PULL  This  program  takes  the  output  supplied  from  HULL  and  provides  the 
mechanism  for  plotting  a  variety  of  graphs,  vector  diagrams  or 
histograms. 

PLANK  This  supplementary  program  is  a  job  monitor.  It  validates  the 

input,  generates  the  various  options  for  the  particular  problem  and 
aborts  if  it  finds  errors  or  inconsistencies. 
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POST  This  supplementary  program  provides  file  management  for  the  whole 
operation.  It  examines  the  input  data  and  then  tailors  the 
composition  of  the  KEEL,  HULL  and  PULL  programs  to  suit. 

In  this  way,  there  is  no  unused  code  incorporated  into  the  run,  thus  ensuring 
that  operating  time  and  cost  are  reduced  to  a  minimum. 


Before  being  able  to  make  an  assessment  of  the  potential  of  a  code 
of  HULL'S  size  and  sophistication,  on  the  Cyber  7600,  some  modifications  to 
the  software  were  naturally  required.  These  changes  were  successfully 
achieved  using  the  sparse  documentation  available. 


3.  THE  PHYSICAL  PROBLEM 


An  outline  of  the  two-dimensional  cartesian  problem  is  shown  in 
Fig.  1.  A  mass  of  shocked  air  (the  shock  wave  is  non -decaying)  travelling 
along  a  tube  is  about  to  impinge  on  a  ramp  angled  at  27s  to  the  horizontal.  A 
small  opening  at  the  top  allows  the  air  flow  to  proceed  beyond  the  wedge.  T 
represents  the  absolute  temperature,  P  the  pressure,  p  the  density  and  u  the 
velocity.  All  symbols  are  stated  in  the  standard  form  according  to  Ben-Dor 
[21. 


The  aim  of  this  work  is  the  generation  of  expertise  with  the  HULL 
code  while  considering  a  problem  that  can  be  examined  in  the  shock  tube  at  the 
Australian  National  University,  Canberra.  Both  approaches  should  show  and 
identify  the  Mach  Stem  and  the  Triple  Point  trajectory. 

The  problem  was  defined  in  KEEL  and  the  data  required  (in  cgs  units! 
are  shown  in  Appendix  A.  The  more  important  details  are  discussed  in  the 
following  paragraphs. 

The  air-filled  tube  is  assumed  to  have  perfectly  reflecting  internal 
walls.  The  equation-of-state  assumes  a  fixed  gamma  law  [3]. 


FIGURE  1.  Problem  Set-up 
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Physical  parameters  for  the  unshocked  and  shocked  air  are  given  in 
Table  1.  “n*e  unshocked  data  are  available  from  KEEL,  while  the  supplied 
shocked-air  values  are  assumed  derived  from  experiment  and  the  input  wave 
velocity,  and  entered  via  the  HULL  data  file  as  shown  in  Appendix  B.  These 
data  will  be  verified  later.  All  other  parameters  are  evaluated  within  HULL. 


TABLE  1. 

PHYSICAL  PARAMETERS  OP  SYSTEM 


QUANTITY 

UNSHOCKED 

AIR 

SHOCKED 

AIR 

Shock  wave  velocity  (m/s) 

- 

5.00877 

X  102 

Ratio  of  specific 

heats 

1 . 4003 

1.4003 

Density  (kg/m3) 

1.225 

2.222 

Particle  velocity 

(m/s) 

0 

2.247  x 

102 

Specific  Internal 

energy  (J/kg) 

2.044  x 

105 

2.692  x 

105 

Pressure  ( Pa ) 

1.013  x 

105 

2.393  x 

105 

Within  KEEL,  a  series  of  monitoring  points  is  set  up  as  shown 
schematically  in  Fig.  2.  These  selected  locations  are  used  to  collect  local 
time-varying  information  of  interest  for  possible  plotting.  Justification  for 
their  use  is  described  later. 


FIGURE  2.  Schematic  showing  station  locations. 

Appendix  C  is  an  example  of  the  input  data  required  for  plotting 
velocity  vectors  and  overpressure  as  a  function  of  time,  at  these  monitoring 


t- 
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stations,  for  the  PULL  program.  Results  for  vector  velocity,  overpressure 
variation  and  pressure  contours  are  given  at  various  times  in  Figures  3,4  and 
5.  The  times  given  at  the  bottom  of  each  plot  correspond  to  a  time  zero  when 
the  shock  front  is  well  before  the  beginning  of  the  ramp. 


In  Figures  3(i)  s  (ii)  the  vector  velocity  of  the  flow  is  shown. 

Here  the  influence  of  the  wedge  on  the  direction  of  flow  is  apparent.  Figs. 
3(iii)  &  (iv)  give  the  velocity  vectors  for  these  same  times,  but  in  a 
logarithmic  scale.  This  can  be  a  useful  display  since  large  changes  to  the 
flow  direction  become  evident  as  can  be  seen,  for  example,  in  Fig.  3(iv)  where 
the  skew  vectors  near  the  front  of  the  flow  by  the  top  boundary  seem  to  give 
early  indication  of  the  influence  of  the  reflective  top  boundary. 


Figures  4(i)  -  4(iv)  show  the  variation  of  pressure  as  a  function  of 
time  at  the  four  locations  A,B,C  and  D  specificed  in  Fig.  2.  These  points 
were  selected  from  theoretical  investigations,  which  indicated  that  the  triple 
point  trajectory  associated  with  a  Mach  Stem  should  pass  between  B  and  C. 
Therefore  measurement  of  a  parameter  such  as  overpressure  should  provide  some 
indication  of  the  program's  applicability. 


If  a  Mach  Stem  were  formed,  as  suggested,  the  overpressure  plots  for 
A  and  B  would  be  significantly  different  from  that  at  C  as  is  observed.  The 
peak  Pm  as  shown  in  Figs.  4(i)  and  4(ii)  can  be  attributed  to  a  Mach  Stem.  In 
Figs.  4(iii)  &  4(iv)  the  first  peak,  P^,  is  caused  by  the  passage  of  the 
incident  shock.  The  second  peak,  Pr,  in  Fig.  4<iii)  &  4(iv)  is  caused  by  the 
shock  reflected  from  the  surface  of  the  wedge  as  it  passes  through  the  pre¬ 
shocked  region  around  the  respective  points  C  &  D  on  its  way  towards  the  top 
boundary.  The  time  interval  between  the  arrival  times  of  the  peaks  P^  and  Pr 
for  locations  C  and  D  and  the  change  in  magnitude  for  Pr  can  be  attributed  to 
the  longer  path  length,  and  therefore  longer  travelling  time  and  greater 
divergence  experienced  by  the  reflected  wave  arriving  at  location  D. 

All  subsequent  peaks  in  Figs.  4(i)  -  4(iv)  are  caused  by  repeated 
reflections  between  the  wedge  and  top  boundary  and  are  of  no  relevance  to  this 
discussion.  Table  2  lists  the  values  of  the  overpressures  P^,  Pr  and  Pffl  at 
all  four  locations. 


Figs.  5(i)-(iii)  show  typical  contour  pressure  plots  obtained  at  the 
indicated  times.  It  would  appear  from  these  plots  that  a  classical  Mach 
Reflection  has  been  found.  Die  incident  and  reflected  waves  are  well  defined 
as  is  the  formation  of  the  Mach  Stem  with  associated  triple  point.  Using 
these  plots  the  position  of  the  triple  point  was  estimated  and  the  trajectory 
found  to  have  an  angle  of  8.2*  ±  0.8°. 
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FIGURE  3(i)  Vector  Velocity  vs  Time 
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FIGURE  3(ii ) 
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FIGURE  4(1)  Overpressure  variation  at  location  A. 
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FIGURE  5(1)  Pressure  contours  showing  Mach  Stem  formation 
along  the  wedge. 
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TABLE  2 


S ELECTED  PEAK  OVERPRESSURE  VALUES  (Pa) 


Locations 

Peak 

P. 

1 

P 

r 

P 

m 

•  A 

- 

- 

2.44  x 

105 

B 

- 

- 

2.45  x 

in 

o 

C 

2.14  x  105 

1.34  x  105 

- 

D 

1.97  x  105 

1.39  x  105 

- 

4.  THEORETICAL  CONSIDERATIONS 


While  there  seems  little  doubt  that  a  Mach  Stem  and  Triple  Point  has 
been  formed  some  calculations  undertaken  during  this  work  were  used  to  monitor 
the  computer  output  and  verify  this  assumption. 

The  Mach  number  of  the  incident  shock,  M^,  is  given  by 
u . 

M.  »  (1) 

i  a 

-  1.472 


where  a  is  the  sonic  velocity  (3.403  x  10^  m/s),  and  u^  is  the  shock  velocity. 

The  associated  pressure  P,  particle  velocity  u,  density  p  and 
specific  internal  energy  E  can  be  derived  from  [4], 


2rM. 


(Y-1) 


Y  +  1 


2(Mi  -  1)a 
(Y  +  DM, 


(2) 
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P  T 

0/0 
o  o 


E 

o 


T 

o 


where  y  is  the  ratio  of  the  specific  heats  and  T  is  the  absolute  temperature 
given  by: 

T,  P,  (Y-1)  (P./P  )  +  (Y+1) 

_L  «  _1  [ - 1 2 - 1 

T  P  (Y+1)  (P./P  )  +  (Y-1) 

O  O  1  o 

and  the  subscripts  o  and  1  refer  to  the  standard  states  (2]  specified  in 
Diagram  1.  Comparison  of  values  calculated  from  (2)  and  those  supplied  to 
HULL  are  given  in  Table  3.  In  each  case,  agreement  is  satisfactory. 


TABLE  3 

VALUES  OF  PARAMETERS  FOR  SHOCKED  AIR 


PARAMETERS 

FROM 

HULL 

INPUT 

FROM 

FORMULAE 

Pressure 

(Pa) 

2.  39  3  X 

105 

2.39  3  X  105 

Particle 

velocity 

(m/s) 

2.247  x 

102 

2.248  x  102 

Density 

(kg/m3) 

2.222 

2.222 

Specific 

Internal 

energy  (J/kg) 

2.692  x 

105 

2.661  X  105 

Now  consider  the  situation  where  the  shock  wave  has  travelled  some 
distance  up  the  inclined  ramp  and  a  Mach  Stem  has  been  formed,  as  indicated  in 
Fig.  6.  Here  is  the  incident  shock  Mach  Number,  Mm  is  the  Mach  Stem  Mach 
Number,  Mr  is  the  reflected  shock  Mach  timber,  9^  the  ramp  angle  and  x  the 
Triple  Point  trajectory  angle. 
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FIGURE  6.  Motion  of  Shock  Front  Up  Ramp. 


Again  using  standard  notation  the  angle  of  incidence  of  the  shock 
wave  with  the  ramp  is  given  by: 

90-0 

w 

whence  to  a  first  approximation,  the  Mach  number  associated  with  the  Mach  Stem 
is  given  by: 


m  sin (90  -  0  ) 
w 

*  1.652 

From  (2)  an  approximation  to  the  pressure  behind  the  Mach  Stem  is: 


2ym  -  (y-1 ) 

_  in  _ 

p  ~  - - - p 

3  T+1  o 


=  3.057  x  10  Pa 


which  leads  to  an  overpressure  value  of  2.044  x  103  Pa.  To  a  first  order, 
this  value  agrees  well  with  the  overpressures  at  Pffl  in  Figs.  4(i)  and  (ii). 

A  more  sophisticated  approach  involves  the  Whitham  theory  [5,6].  By 
using  an  auxiliary  function  A  defined  by 

A  -  f  (M) 


s  exp  {-{in  ^  In  (M2  -  +  In  j-. 


*  1/2  tu  +  1/2 1  ■  1/2  in  tn  -  <^jr> 1/2 1 


+  ( 


V2 


Y(Y-1) 


)  '*  In  [(M2  +  tA->  1/2  +  (M2  -  1/2 


Y-1 


2Y 


*  'ttW'  1/2 '“’’i 


■r  -  ,,^,l2l  mV:  4  r~" — <*> 

4Y  /2  (Y-1)  IM  +  <  2/  (Y-1 )  )  1  /2  [M  -  (Y-D/2Y1  /2 


where  u 


2  _  (Y-1)  M  +  2 


(5) 


2Y  M  -  (Y-1) 

values  of  the  ramp  angle  and  the  triple  point  trajectory  can  be  derived  from 


tan  0 

w 


tan  x 


(M 


M.2)  1/2 


(A.  -  A 

l  m 


y/2 


A  M  +  A.M. 
mm  li 


A  1 

—  (  - 

A  1  2 

1  1  -  (A  /A. r 


(M./M  )  v 
1  hi  j  72 


m 


(Pv 


(7) 


Using  the  initial  estimate  for  Mm  of  1.652  and  the  Known  value  of 
Mi#  values  of  Am  and  Ai  can  be  calculated  from  (4)  and  (5),  and  hence  0W  and 
X  can  be  calculated  from  (6)  and  (7).  In  fact,  if  ©w  is  assumed  constant,  a 
better  value  for  can  be  obtained  and  the  procedure  repeated.  On  successful 
completion  of  this  iterative  process,  (2)  can  be  used  to  calculate  the  Mach 
Stem  pressure  and  thus  the  overpressure.  Results  are  shown  in  Table  4. 


TABLE  4 


SHOCK  PARAMETERS  FROM 

WHITHAM 

THEORY 

Mach  Number  of  Mach  Stem 

1 .778 

Angle  of  Ramp 

0w 

27° 

Triple  Point  Trajectory 

X 

8° 

Overpressure  behind  Mach  Stem 

P 

2.56  : 

The  overpressure  behind  the  Mach  Stem  as  derived  and  shown  in  Table 
4  is  in  very  good  agreement  with  that  computed  from  locations  A  and  B  and 
shown  in  Table  2. 
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The  Whitham  theory  can  also  be  used  to  predict  the  Mach  number  Mr  of 
the  wave  as  reflected  at  the  ramp  given  by 


((Y  +  1/2Y)  [M  -  (Y-1  )/2y  ] 

XzJ  .  _ 52 _ 

2y  +  2 

M.  -  (y  -1  )/2y 


1.421 


or  M  =  1.192 

r 


The  velocity  of  this  reflected  wave  is  not  constant;  its  overpressure  varies 
depending  on  both  the  path  length  and  instantaneous  velocity.  A  range  for  the 
expected  overpressures  is  between  that  corresponding  to  the  Mach  number  of 
1.192  (hence  a  maximum  overpressure  of  2.56  x  105  Pa)  and  that  of  the  incident 
shock  of  1.38  x  10s  Pa.  Comparison  with  the  HULL  values  for  peak  Pr  at  C  and 
D  on  Table  2  confirms  this  result. 


Finally,  the  theoretical  value  for  the  Triple  Point  trajectory  (as 
shown  in  Table  4)  is  in  good  agreement  with  the  value  of  8.2°  as  measured  from 
Figures  5(i)  -  5(iii). 


5.  CONCLUSION 


The  sophisticated  computer  code  HULL  can  be  used  to  simulate  a  shock 
front  impinging  on  an  inclined  ramp.  Many  more  complex  investigations  are 
planned  both  within  a  shock  tube  and  in  the  study  of  blast  phenomena  from 
detonating  explosive  charges.  The  next  phase  in  this  work  will  incorporate 
the  development  of  a  User  Manual  for  the  Australian  version  of  the  code. 
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comments  during  the  preparation  of  this  Report. 
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APPENDIX  A 


Input  Data  for  KEEL  run 


♦EOS 

KEEL 

PROB  4.0 

BBOUND=0 

TBOUND-O 

LB0UND=7 

RBOUND-1 

ATMOS=5 

COLD=TRUE 

DIMEN-2 

EOS-2 

GEOM-1 

HEADER 

STEP  SHOCK  UP  A  27  DEGREE  RAMP 
ISLAND-1 
SHORE=1 
IMAX-300 
JMAX-72 
T=*0.0 

PTSTOP-0.080 

REZONE-O 

NSTN-4 

RELGAM-1 .4003 
RELRHO-1 . 225E-3 
RELSIE-2.044E9 
RELPO-1.01325E6 
MESH 

XO-O.O  NX-300  DX-0.25 

YO— 0.0  NY-  72  DY-0.25 

GENERATE  PACKAGE  AIR 

I-0.269152E10 

U-0.224738E5 

V-0.0 

RHO-0. 002222 
RECTANGLE 
XI— 3000  X2-5.0 

Y1— 3000  Y2-3000 

PACKAGE  SHORE  TRIANGLE 

X1-11.00  Y1-0.0  X2-45.0  Y2-17.0  X3-45.0  Y3-0.0 

PACKAGE  ISLAND 

RECTANGLE 

X1-45.0  X2-75.0 

Y1-0.0  Y2-17.0 

STATIONS 

XS-21.0  YS-11.0 

XS-33.0  YS-12.0 

XS— 33.0  YS— 14.0 
XS-33.0  YS-16.0 

♦EOP 


APPENDIX  B 


Input  Data  for  HULL  run 


*EOS 

HULL 

PROB  4.0000 
CYCLE-440 
INPUT 
T-6.0E-4 
SHKTST— 9.982E-5 
SHKAZM-0. 0 
SHKWVL— 50087. 7 
SHKPVL— 22473.8 
SHKEBS— 2. 691 52E9 
SHKRHO-O. 002222 
MRELER— 1 . OE-6 
REZONE— 0 
STABF-0.5 
CST0P-700 
*BOP 


APPENDIX  C 


Input  Data  for  PULL  Run  to  Plot: 


(A)  Vector  Velocity 


*EOS 

PULL 

PROB  4.0 

PL0TPKG=6 

PLTXPG=7 

PLTYPG=10 

PLTFTR=5 

WECT 

CTIME=0. 001 ,0.0011,0.0012 
*E0P 


(B)  Station  Data 


*E0S 

PULL 

PROB  4.0 

PLOTPKG-6 

PLTXPG-7 

PLTYPG-10 

STATION 

PLTFTRat  1 0 

PSTA-19 

LSTA*24 

*EOP 
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